Intensity and shapefile data is handled here.
# make dates
info = gdalinfo("VV_2017_clip.tif")
descr = info[grepl(info, pattern = "Description = ")]
descr = gsub(descr, pattern = "Description = ", replacement = "")
dates2017 <- as_date(descr)
info = gdalinfo("VV_clip.tif")
descr = info[grepl(info, pattern = "Description = ")]
descr = gsub(descr, pattern = "Description = ", replacement = "")
dates2018 <- as_date(descr)
# function takes a shape as an AOI and a year (2017 or 2018)
loadPolygon <- function(shape, year) {
if(year == 2017) {
int_VV <- int_VV_2017
int_VH <- int_VH_2017
dates <- dates2017
} else {
if(year == 2018) {
int_VV <- int_VV_2018
int_VH <- int_VH_2018
dates <- dates2018
} else {
return("wrong year entered!")
}
}
int_vv <- st_as_stars(int_VV[shape])
int_vh <- st_as_stars(int_VH[shape])
int_vv <- st_set_dimensions(int_vv, 3, val = dates, names = "time")
int_vh <- st_set_dimensions(int_vh, 3, val = dates, names = "time")
# c()
comb <- c(int_vv, int_vh, along="bands")
# switch bands to attributes, this is more of a design choice.
comb_split <- split(comb, "bands")
names(comb_split) <- c("VV", "VH")
return(comb_split)
}
area1 <- loadPolygon(shape[124,], 2018)
area2 <- loadPolygon(shape[124,], 2017)
# this kind of water puddle:
plot(area1[1,,,1])
Before, thresholds were calculates for single years, resulting in -15.3 for 2018 and about -18 for 2017.
# prepare data
val.svm$type <- as.factor(val.svm$type)
attach(val.svm)
x <- subset(val.svm, select=c(-type, -sitch))
y <- type
# make model
svmmod <- svm(x,y)
# summary(svmmod)
# make prediction, confusion matrix
pred <- predict(svmmod,x)
table(pred, y)
## y
## pred field water
## field 4019 68
## water 1 1948
# threshold is
svmmod$x.scale$`scaled:center`
## [1] -17.62537
# plot classified with threshold calculated by SVM
plot(area1[1,,,1], col = c("white", "black"), breaks = c(-35, svmmod$x.scale$`scaled:center`, 10))
# vs old threshold
# plot(all[1,,,1], col = c("white", "black"), breaks = c(-35, thresh$threshold, 10))
ggplot(val.svm, aes(x=sitch, y=VV)) +
geom_boxplot() +
geom_hline(yintercept = svmmod$x.scale$`scaled:center`)
## The following objects are masked from val.svm (pos = 3):
##
## sitch, type, VV
## y
## pred field water
## field 4060 130
## water 5 2000
## [1] -17.55717
allRas2017 <- raster("./minus17/allRas2017.tif")
allRas2018 <- raster("./minus17/allRas2018.tif")
colo = viridisLite::inferno(29)
breaks = seq(0, 29, 1)
image(allRas2017, col = colo, breaks = breaks, main="2017")
colo = viridisLite::inferno(31)
breaks = seq(0, 31, 1)
image(allRas2018, col = colo, breaks = breaks, main="2018")
allRas2017 <- raster("./allRas2017.tif")
allRas2018 <- raster("./allRas2018.tif")
colo = viridisLite::inferno(29)
breaks = seq(0, 29, 1)
image(allRas2017, col = colo, breaks = breaks, main="2017")
colo = viridisLite::inferno(31)
breaks = seq(0, 31, 1)
image(allRas2018, col = colo, breaks = breaks, main="2018")
allRas2017 <- read_stars("./minus17/allRas2017.tif")
allRas2018 <- read_stars("./minus17/allRas2018.tif")
st_crs(allRas2017) <- "EPSG:25832"
st_crs(allRas2018) <- "EPSG:25832"
colo = viridisLite::inferno(29)
breaks = seq(0, 29, 1)
image(allRas2017[study_area[1,]], col = colo, breaks = breaks, main="2017")
colo = viridisLite::inferno(31)
breaks = seq(0, 31, 1)
image(allRas2018[study_area[1,]], col = colo, breaks = breaks, main="2018")
st_crs(allRas2018) <- "EPSG:25832"
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
# cut shapefile to area of interest
inter_shape <- shape[lengths(st_intersects(shape, study_area[1,])) != 0,]
plot(inter_shape)
# load
moor2017 <- loadPolygon(inter_shape, 2017)
plot(moor2017)
# set - 23 as Threshold for open water
moor2017[moor2017[1,,,] < -23.5] <- NA
# plot(moor2017)
makeTimeSeries <- function(stars, shape) {
# for each shape
for (i in 1:length(shape$Förderung)) {
cut <- stars[shape[i,]]
cut <- as.data.frame(st_apply(cut, "time", mean, na.rm = TRUE))
if(!is.na(shape[i,]$Förderung)) {
names(cut) <- c("time", paste0("F", i, "VV"), paste0("F", i, "VH"))
} else {
names(cut) <- c("time", paste0(i, "VV"), paste0(i, "VH"))
}
# at first iteration
if(i == 1) {
# create master df
df <- cut
# at every ther iteration
} else {
# append
df <- cbind(df, cut[,2:3])
}
}
return(df)
}
plotTimeSeries <- function(df, year) {
if(year==2017) {
rain.select <- rain.all.2017.df
} else {
if(year==2018) {
rain.select <- rain.all.2018.df
}
else {
if(class(year) == "data.frame") {
rain.select <- year
}
}
}
g <- ggplot() + geom_bar(aes(x=rain.select$time, y=rain.select$mean), stat='identity', fill = "red", alpha = 0.8) + coord_cartesian(ylim=c(0,35)) + ggtitle("Mean of Moorflächen (*2) plotted against Precipitation (3 day sum) \nblue = in Förderung, black = nicht in Förderung") + xlab("Time") + ylab("Precipitation in mm, Backscatter VV")
#+ scale_y_continuous(sec.axis = sec_axis(~. /10 + 22, name = "Intensity + 22"))
# select VV columns
for (i in seq(2, ncol(df), 2)) {
if(names(df)[i] == paste0("F", i/2, "VV")) {
g <- g + geom_line(aes_string(x=df[,1], y = df[,i]*2 + 44, alpha = 0.7), color = "blue")
g <- g + geom_point(aes_string(x=df[,1], y = df[,i]*2 + 44, alpha = 0.7), color = "blue")
} else {
g <- g + geom_line(aes_string(x=df[,1], y = df[,i]*2 + 44), alpha = 0.3)
}
}
print(g)
}
df <- makeTimeSeries(moor2017, inter_shape)
plotTimeSeries(df, 2017)
# load
moor2018 <- loadPolygon(inter_shape, 2018)
# set - 23 as Threshold for open water
moor2018[moor2018[1,,,] < -23.5] <- NA
plotTimeSeries(makeTimeSeries(moor2018, inter_shape), 2018)
dff <- rbind(makeTimeSeries(moor2017, inter_shape), makeTimeSeries(moor2018, inter_shape))
plotTimeSeries(dff, year = rbind(rain.all.2017.df, rain.all.2018.df))
dfff <- dff[,c(1,seq(2,ncol(dff), 2))]
gef <- c(1)
ngef <- c(1)
for (i in 2:length(names(dfff))) {
str <- paste0("F", i - 1, "VV")
if(names(dfff)[i] == str) {
gef <- c(gef, i)
} else {
ngef <- c(ngef, i)
}
}
inF <- dfff[,gef]
ninF <- dfff[,ngef]
inF$Mean <- rowMeans(inF[,2:ncol(inF)])
ninF$Mean <- rowMeans(ninF[,2:ncol(ninF)])
dfM <- cbind(inF[,c(1, ncol(inF))], ninF$Mean, ninF$Mean)
names(dfM) <- c("time", "F1VV", "1XX", "1VV")
plotTimeSeries(dfM, year = rbind(rain.all.2017.df, rain.all.2018.df))
## Warning in if (year == 2017) {: the condition has length > 1 and only the first
## element will be used
## Warning in if (year == 2018) {: the condition has length > 1 and only the first
## element will be used
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
head(dfM)
## time F1VV 1XX 1VV
## 1 2017-01-08 -14.97066 -15.63378 -15.63378
## 2 2017-02-01 -12.78420 -12.97320 -12.97320
## 3 2017-02-13 -15.60427 -14.87243 -14.87243
## 4 2017-02-25 -11.78134 -12.13760 -12.13760
## 5 2017-03-09 -11.77774 -12.27017 -12.27017
## 6 2017-03-21 -11.97166 -11.94932 -11.94932
# delete the column needed for plotting function
dfM <- dfM[,c(1,2,4)]
# all acquisitions are 12 days apart, the first one is 8.1.2017
# rain
rain.all <- rbind(rain.all.2017.df, rain.all.2018.df)
# make sum rain vector
vec <- c(rain.all[1,2] + rain.all[2,2] + rain.all[3,2] + rain.all[4,2] + rain.all[5,2] + rain.all[6,2] + rain.all[7,2] + rain.all[8,2])
for (i in 1:nrow(dfM)-1) {
offset <- 8
inter <- 12
count <- i - 1
start <- count*inter + offset + 1
end <- i * inter + offset
sum <- 0
# print(c(start, end))
for (j in start:end) {
sum <- sum + rain.all[j,2]
}
vec <- c(vec, sum)
}
df <- cbind(dfM, vec)
names(df) <- c("time", "FVV", "nFVV", "prec")
ggplot(df, aes(x = time)) +
geom_bar(aes(x = time - 6, y = prec), stat = 'identity', fill = "lightblue", alpha = 0.8) +
geom_line(aes(y = FVV * 4 + 70, color = "in Förderung")) +
geom_point(aes(y = FVV * 4 + 70, color = "in Förderung")) +
geom_line(aes(y = nFVV * 4 + 70, color = "nicht in Förderung")) +
coord_cartesian(ylim = c(0,60)) +
ggtitle("Mean of Moorflächen and Precipitation") + xlab("Time") +
ylab("Precipitation in mm/m²") +
scale_y_continuous(sec.axis = sec_axis(~. *0.25 -17.5, name = "Intensity in dB")) +
scale_color_manual(name = "Förderung", values = c("in Förderung"="blue", "nicht in Förderung"="black")) +
theme(legend.position = "bottom")
# make sum rain vector
vec <- c(rain.all[1,2] + rain.all[2,2] + rain.all[3,2] + rain.all[4,2] + rain.all[5,2] + rain.all[6,2] + rain.all[7,2] + rain.all[8,2])
for (i in 1:nrow(dfM)-1) {
offset <- 8
inter <- 12
count <- i - 1
start <- count*inter + offset + 1
end <- i * inter + offset
sum <- 0
# print(c(start, end))
weight <- 1
for (j in start:end) {
sum <- sum + (1/weight * rain.all[j,2])
weight <- weight + 1
}
vec <- c(vec, sum)
}
gf <- cbind(dfM, vec)
names(gf) <- c("time", "FVV", "nFVV", "prec")
ggplot(gf, aes(x = time)) +
geom_bar(aes(x = time - 6, y = prec), stat = 'identity', fill = "lightblue", alpha = 0.8) +
geom_line(aes(y = FVV * 1.5 + 28, color = "in Förderung")) +
geom_point(aes(y = FVV * 1.5 + 28, color = "in Förderung")) +
geom_line(aes(y = nFVV * 1.5 + 28, color = "nicht in Förderung")) +
coord_cartesian(ylim = c(0,25)) +
ggtitle("Mean of Moorflächen and Precipitation, Weight: 1 / # of Day") + xlab("Time") +
ylab("Precipitation in mm/m²") +
scale_y_continuous(sec.axis = sec_axis(~. *2/3 - 56/3, name = "Intensity in dB")) +
scale_color_manual(name = "Förderung", values = c("in Förderung"="blue", "nicht in Förderung"="black")) +
theme(legend.position = "bottom")
head(df)
## time FVV nFVV prec
## 1 2017-01-08 -14.97066 -15.63378 15.3018018
## 2 2017-02-01 -12.78420 -12.97320 10.7297298
## 3 2017-02-13 -15.60427 -14.87243 5.3207207
## 4 2017-02-25 -11.78134 -12.13760 0.8630631
## 5 2017-03-09 -11.77774 -12.27017 22.7909910
## 6 2017-03-21 -11.97166 -11.94932 22.4567567
dff <- df
dff$FVVbyprec <- (25 + dff$FVV - dff$prec) / (25 + dff$FVV + dff$prec)
ggplot(dff, aes(x = time)) +
geom_line(aes(x = time - 6, y = FVVbyprec)) +
geom_bar(aes(x =time - 6, y = prec), stat = 'identity', fill = "lightblue", alpha = 0.8) +
geom_line(aes(y = FVV * 4 + 70, color = "in Förderung")) +
geom_point(aes(y = FVV * 4 + 70, color = "in Förderung")) +
geom_line(aes(y = nFVV * 4 + 70, color = "nicht in Förderung")) +
coord_cartesian(ylim = c(-5,40)) +
ggtitle("Mean of Moorflächen and Precipitation") + xlab("Time") +
ylab("Precipitation in mm/m²") +
scale_y_continuous(sec.axis = sec_axis(~. *0.25 -17.5, name = "Intensity in dB")) +
scale_color_manual(name = "Förderung", values = c("in Förderung"="blue", "nicht in Förderung"="black")) +
theme(legend.position = "bottom")
dff$Fchange <- 1:nrow(dff)
for (i in 2:nrow(dff)) {
dff[i,6] <- (dff[i,2] - dff[i-1,2]) / (dff[i,2] + dff[i-1,2])
}
dff$Rchange <- 1:nrow(dff)
for (i in 2:nrow(dff)) {
dff[i,7] <- (dff[i,4] - dff[i-1,4]) / (dff[i,4] + dff[i-1,4])
}
dff$FRchange <- 1:nrow(dff)
for (i in 2:nrow(dff)) {
dff[i,8] <- (dff[i,6] - dff[i-1,7]) / (dff[i,6] + dff[i-1,7])
}
# normalized change of normalized change over time with previous time step
ggplot(dff, aes(x = time)) +
geom_line(aes(x = time - 6, y = FRchange)) +
geom_bar(aes(x =time - 6, y = prec), stat = 'identity', fill = "lightblue", alpha = 0.8) +
geom_line(aes(y = FVV * 4 + 70, color = "in Förderung")) +
geom_point(aes(y = FVV * 4 + 70, color = "in Förderung")) +
geom_line(aes(y = nFVV * 4 + 70, color = "nicht in Förderung")) +
coord_cartesian(ylim = c(-5,40)) +
ggtitle("Mean of Moorflächen and Precipitation") + xlab("Time") +
ylab("Precipitation in mm/m²") +
scale_y_continuous(sec.axis = sec_axis(~. *0.25 -17.5, name = "Intensity in dB")) +
scale_color_manual(name = "Förderung", values = c("in Förderung"="blue", "nicht in Förderung"="black")) +
theme(legend.position = "bottom")
# make cumsum
dff$Rsum <- cumsum(dff$prec - mean(dff$prec[1:30]))
# start from mean prec
dff$prec_1 <- c(mean(dff$prec[1:30]), dff$prec[2:60])
dff$Rsum_1 <- cumsum(dff$prec_1 - mean(dff$prec[1:30]))
# calc means
mean(dff$prec) # all
## [1] 19.53345
mean(dff$prec[1:30]) # 2017
## [1] 26.89676
mean(dff$prec[31:60]) # 2018
## [1] 12.17015
dff$Rsum <- dff$Rsum - mean(dff$prec)
# ggplot(dff[1:30,], aes(x = time)) +
# geom_bar(aes(y = Rsum / 50), stat = 'identity') +
# geom_line(aes(y = FVV + 11.5)) +
# ggtitle("cumsum rain, scaled")
# ggplot(dff[1:30,], aes(x = time)) +
# geom_line(aes(y = Rsum / 50), color="lightblue") +
# geom_line(aes(y = FVV + 11.5)) +
# ggtitle("cumsum rain, scaled")
# R sum 1
ggplot(dff[1:30,], aes(x = time)) +
geom_bar(aes(y = Rsum_1 / 50), stat = 'identity') +
geom_line(aes(y = FVV + 12)) +
ggtitle("cumsum rain, scaled, starting at 1 = mean")
ggplot(dff[1:30,], aes(x = time)) +
geom_line(aes(y = Rsum_1 / 50), color="lightblue") +
geom_line(aes(y = FVV + 12)) +
ggtitle("cumsum rain, scaled, starting at 1 = mean")
# weighted stuff
gf$Rsum <- cumsum(gf$prec - mean(gf$prec[1:30]))
ggplot(gf[1:30,], aes(x = time)) +
geom_bar(aes(y = Rsum / 10), stat = 'identity') +
geom_line(aes(y = FVV + 11.5)) +
ggtitle("cumsum weighted rain, 2017")
ggplot(gf[1:30,], aes(x = time)) +
geom_line(aes(y = Rsum / 10), color="lightblue") +
geom_line(aes(y = FVV + 11.5)) +
ggtitle("cumsum weighted rain, 2017")
# 2018
dff$Rsum[31:60] <- cumsum(dff$prec[31:60] - mean(dff$prec[31:60]))
ggplot(dff[31:60,], aes(x = time)) +
geom_line(aes(y = Rsum / 4), color="lightblue") +
geom_line(aes(y = FVV + 17)) +
ggtitle("cumsum rain 2018, scaled, - mean 2018")
ggplot(dff[1:30,], aes(x = time)) +
geom_bar(aes(x =time - 6, y = prec), stat = 'identity', fill = "lightblue", alpha = 0.8) +
geom_line(aes(y = FVV * 4 + 70, color = "in Förderung")) +
geom_point(aes(y = FVV * 4 + 70, color = "in Förderung")) +
geom_line(aes(y = nFVV * 4 + 70, color = "nicht in Förderung")) +
geom_hline(yintercept = mean(dff$prec[1:30]), color = "lightblue") +
coord_cartesian(ylim = c(-5,40)) +
ggtitle("Mean of Moorflächen and Precipitation") + xlab("Time") +
ylab("Precipitation in mm/m²") +
scale_y_continuous(sec.axis = sec_axis(~. *0.25 -17.5, name = "Intensity in dB")) +
scale_color_manual(name = "Legend", values = c("in Förderung"="blue", "nicht in Förderung"="black", "Mean of Precipitation 2017"="lightblue")) +
theme(legend.position = "bottom")
# dff$prec_2 <- c(mean(dff$prec), dff$prec[2:60])
dff$Rsum_2 <- cumsum(dff$prec - 23)
ggplot(dff, aes(x = time)) +
geom_bar(aes(y = Rsum_2 / 50), stat='identity') +
geom_line(aes(y = FVV * 1.2 + 15))
acf(dff$prec)
pacf(dff$prec)
acf(dff$FVV)
pacf(dff$FVV)
cor(dff$FVV, dff$prec)
## [1] 0.256756
ccf(dff$FVV, dff$prec)
ccf(dff$FVV, dff$Rsum)
ccf(dff$FVV, dff$Rsum_1)
acf(dff$Rsum_1)
dff$acq <- 1:nrow(dff)
f <- function(x) sum((dff$FVV - (x[1] + x[2] * sin(pi * (dff$acq+x[3])/30)))^2)
nlm(f,c(0,0,0))
## $minimum
## [1] 82.44628
##
## $estimate
## [1] -13.386189 1.610426 -8.671524
##
## $gradient
## [1] -1.111501e-06 -1.299817e-05 7.397523e-06
##
## $code
## [1] 1
##
## $iterations
## [1] 12
dff$FVV.per <- -13.386 + 1.61 * sin(pi*(dff$acq - 8.6715)/30)
ggplot(dff, aes(x = time)) +
geom_line(aes(y = FVV)) +
geom_line(aes(y = FVV.per, color = "red"))
ggplot(dff, aes(x = time)) +
geom_bar(aes(y = prec / 10), stat = 'identity') +
geom_line(aes(y = FVV - FVV.per + 1.8))
an <- dff$FVV - dff$FVV.per
an.ar5 = arima(an, c(5,0,0))
an.ar5
##
## Call:
## arima(x = an, order = c(5, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 intercept
## 0.3567 0.2003 -0.2136 -0.0615 -0.0723 0.0182
## s.e. 0.1295 0.1439 0.1624 0.1626 0.1580 0.1705
##
## sigma^2 estimated as 1.056: log likelihood = -87, aic = 187.99
acf(an, type="partial")
arima(an, c(1,0,0))$aic
## [1] 185.1914
dfh <- dff[1:30,]
n <- 4
f <- function(x) sum((dfh$FVV - (x[1] + x[2] * sin(pi * (dfh$acq+x[3])/n)))^2)
mod <- nlm(f,c(0,0,0))
dfh$FVV.per <- mod$estimate[1] + mod$estimate[2] * sin(pi*(dfh$acq + mod$estimate[3])/n)
ggplot(dfh, aes(x = time)) +
geom_line(aes(y = FVV)) +
geom_line(aes(y = FVV.per, color = "red"))
n <- 30
f <- function(x) sum((dff$prec - (x[1] + x[2] * sin(pi * (dff$acq+x[3])/n)))^2)
mod <- nlm(f,c(0,0,0))
dff$prec.per <- mod$estimate[1] + mod$estimate[2] * sin(pi*(dff$acq + mod$estimate[3])/n)
ggplot(dff, aes(x = time)) +
geom_line(aes(y = FVV)) +
geom_line(aes(y = prec.per, color = "prec")) +
geom_line(aes(y = FVV.per * 10 + 150, color="FVV")) +
ggtitle("Annual Trends, scaled")
ggplot(dff, aes(x = time)) +
# geom_bar(aes(y = prec / 10, color = "rain"), stat = 'identity') +
geom_bar(aes(y = (prec - prec.per) / 10, color = "residuals"), stat = 'identity') +
geom_line(aes(y = FVV - FVV.per + 1.8)) +
ggtitle("residual rain and residual VV")